Comparing the first and second order tlieories of relativistic dissipative fluid dynamics 
using the 1+1 dimensional relativistic flux corrected transport algorithm 
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Focusing on the numerical aspects and accuracy we study a class of bulk viscosity driven expansion 
scenarios using the relativistic Navier-Stokes and truncated Israel-Stewart form of the equations of 
relativistic dissipative fluids in 1+1 dimensions. The numerical calculations of conservation and 
transport equations are performed using the numerical framework of flux corrected transport. We 
show that the results of the Israel-Stewart causal fluid dynamics are numerically much more stable 
and smoother than the results of the standard relativistic Navier-Stokes equations. 
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I. INTRODUCTION 

The recent discovery of near perfect fluidity of hot 
QCD matter at the Relativistic Heavy Ion Collider 
(RHIC) 1] brought a lot of attention and interest in mod- 
eling the collective phenomena in relativistic heavy-ion 
collisions using the relativistic dissipative fluid dynamical 
approach. In contrast to perfect fluid dynamical models, 
dissipative fluids provide a more accurate and physically 
more plausible description incorporating first and second 
order corrections compared to perfect fluids. 

These higher order corrections are irreversible; ther- 
mal conductivity and dissipation, related to temperature 
gradients and inhomogeneities of the flow field. A linear 
relation between the two establishes transport equations, 
where the parameters entering these equations are the so- 
called transport coefficients, for thermal conductivity A, 
shear viscosity 77, and bulk viscosity C, also referred to as 
second viscosity or volume viscosity @, y, Q ■ 

Recently many studies have specifically investigated 
the fluid dynamical description of matter created at 
RHIC including shear viscosity, see 0, i, 0, S, i, [M [IH, 

fand references therein. Some of these calculations 
[lOl | made use of the first order theory by Eckart [2], 
and Landau and Lifshitz [3|, but the main focus was on 
the second order causal theory of dissipative fluid dynam- 
ics by Israel and Stewart [4| , and the theory by Ottinger 
and Grmela [l3l |. 

These calculations particularly examined the effect of a 
small shear viscosity motivated by the conjectured lower 
bound from the AdS/CFT correspondence [llllllll. 
It was found only recently that, contrary to perturbative 
QCD estimates [Ij], lattice QCD reveals a large increase 
of the bulk viscosity near the critical temperature [1^, 
IiqI . [201 ■ This numerical evidence motivates studies of the 
evolution of matter with large viscosity. It has also been 
suggested that a large bulk viscosity near Tc may entirely 
change the standard picture of adiabatic hadronization 
employed so far in hydrodynamical models [2l| . 

In this paper we address the phenomena related to 
viscous evolution of matter in 1+1 dimensional systems 
neglecting the contribution of heat conduction. We fo- 



cus on the numerical implementation of both the 1st- 
order and 2nd-order approaches and investigate specific 
test cases to clarify numerical aspects and accuracy of 
the solutions. This represents an important first step 
before multi-dimensional models can be constructed and 
applied. We also show how the relaxation equations for 
the dissipative corrections in the 2nd-order theory can 
be solved efficiently and accurately also via the flux cor- 
rected transport algorithm by writing them in the form 
of continuity equations with a source. 

The paper is organized as follows. First, we briefly 
recapitulate and formulate the equations of dissipative 
fluid dynamics and the numerical method which will be 
used to solve the respective equations. Afterwards we 
present and discuss the results in several cases. To our 
knowledge major parts of this work including the speciflc 
comparisons between the Ist-order and 2nd-order theo- 
ries is presented and discussed in detail for the first time. 



II. THE EQUATIONS OF DISSIPATIVE FLUID 
DYNAMICS 



We adopt the standard notation for four-vectors and 
tensors and use the natural units H = c = k = 1, through 
this paper. The upper greek indices denote contravariant 
while the lower indices denote covariant objects. The ro- 
man indices or the bold face letters denote three- vectors. 

In the Eckart frame, the conserved charge four-current 
is, N^ — nu^, where n is the local rest frame con- 
served charge density, u^ — (7, 7v) is the four-fiow of 
matter normalized to one, u'^u^ = 1, and the relativis- 
tic gamma is, 7 = l/yl^-v^. The dissipative energy- 
momentum tensor is, T^'' ~ {e + p + iVju'^u'^ ^ (p + 
iTjgf^" + tt'"', where, e = u^.T^'^u^, is the local rest 
frame energy density, the orthogonal projection of the 
energy-momentum tensor, p(e, n) + H = — ■iAp^T'*'', de- 
notes the local equilibrium pressure plus the bulk pres- 



sure and g^^" = g^^, — diag(l,— 1 


,— 1,— 1), is the met 


ric of the fiat space-time. The stress tensor, tt'"' = 


j (a^a;^ + aj;a^) - iA'^-A^^ 


T"'^, is the symmet 



ric, traceless ir'^'^g^i, = 0, and orthogonal to the flow ve- 



locity, TT^^Ui, = 0, part of the energy-momentum tensor. 
The local conservation of charge, energy, and momentum 
requires that. 






0, 
0, 



(1) 
(2) 



and the second law of thermodynamics demands that 
the four-divergence of the entropy four-current is non- 
decreasing and positive. 



5^^^ > 0. 



(3) 



Here 9^ = {dt,di) denotes the four-divergence where 
dt = d/dt is the time-derivative and di = d/dxi — 
{d/dx, d/dy, d/dz) is the divergence operator. 

The explicit form of the conservation equations for 
charge, energy, and momentum are. 



^tN''+^^N'' = 0, 


(4) 


dtT°° + d,T''' = 0, 


(5) 


dtT"^ + d,T'^ = 0, 


(6) 



where we defined the conserved charge density, N^, 
charge flux, iV% the total energy density, T"", the en- 
ergy flux density, T"*, the momentum flux density, T*'^ 
and the momentum flux density tensor, T*-' . These labo- 
ratory frame quantities can be expressed in terms of the 
local rest frame quantities and velocity as. 
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717, 

{e+p + U)-/^ -{p + n) 

{e+p + Uyf^v, + TT°\ 

WzT™ + w,(p + n) - i;,^o" + ^OV 

„2„ 



(7) 

(8) 

(9) 

(10) 



T'^ = {e+p + U)Yv^Vj -{p + n).g'J + n'^ , (11) 



f,T°J" - (p + n)g'- 



rOj" 



The relation between the local rest frame and labora- 
tory frame quantities can be calculated using the above 
equations, hence 



n 
e 









(12) 
(13) 



where the absolute value of the velocity is, v = |v|. These 
local rest frame quantities are needed to calculate the 
pressure, p{e,n), from the equation of state (EOS). 

The fluid velocity and relativistic gamma can be cal- 
culated from eq. pTI)l . therefore. 



(TO' - ttO*) 



fJ^OO 

1 



P(e,n,n) 
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vT" 



(14) 
(15) 



where P(e, n, 11) = p{e, n) +11, gives the correction to the 
equilibrium pressure absorbed in the trace of the energy- 
momentum tensor. 



A. 1-|-1 dimensional expansion 

For simple 14-1 dimensional systems in 1+3 dimen- 
sional space-time, where u^ = 7(1, 0, 0, w^), the equations 
of dissipative fluid dynamics reduce to a similar form as 
in the case of perfect fluids. Let us denote the pres- 
sure in the longitudinal direction by, Pz=p + Il + Tr = 
P(e,n,n) + TT, where tt = 7r^^/7^ is the local rest frame 
value of the stress. Due to construction, the traceless- 
ness property implies that tt^^ = tt*'*' = —tt/2, and 
7r°'' — w^7^7r. The orthogonality relations will further re- 
duce the number of unknowns; note that, 7r°^ = 7r°^ — 0, 
^Oz ^ ■y^TT^^ = 11^7271-^ and all non-diagonal components 
of the stress tensor vanish, Tr*;'^ ■ = 0, thus the only com- 
ponent of the shear tensor we have to propagate is n [22] ■ 

The conservation equations follow from eqs. (I4l5l6p . 



= 0, 



9tT°" + a,(t>,T°°) 

where the laboratory frame quantities are, 

A^" = ^7, 

TOO ^ (e + F,)72_ 

The local rest frame variables expressed through the lab- 
oratory frame quantities, the velocity and relativistic 
gamma are. 



dz{v,P,), 


(17) 


dzPz, 


(18) 


;ies are. 






(19) 
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(20) 
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(21) 



n = N^^/l-v 

rpOz 
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while the local rest frame effective pressure is, 

Pz = p{e, n) + n + TT . 



(22) 
(23) 

(24) 
(25) 

(26) 



Here the equilibrium pressure is given by the equation of 
state, p — c^e, where Cs is the local speed of sound. The 
equations and quantities for a perfect fluid are obtained 
in the limit of vanishing dissipation corresponding to, 
Pz — !■ p(e,n), while the form of conservation equations 
and the expressions relating the laboratory frame quan- 
tities to the rest frame quantities and the calculation of 
the velocity are formally the same as for perfect fluids. 

The last two variables that remain to be explicitly de- 
fined are the bulk pressure and the shear. These can be 
calculated from eqs. (|1I2I3P using different approaches. 
To study the various methods is out of the scope of the 
current manuscript, however these theories and recent 
new phenomenological development aimed to extend the 
theory of dissipative fluids sheds light on the open ques- 
tions related to the ambiguities on this matter, see for 



example [23|, |2J, |25|, 123, 1211, 123, |29|, |30|, IM 121 and refer- 
ences therein. 

In the Ist-order theories of Eckart [2] or Landau and 
Lifshitz 3], i.e., the relativistic Navier-Stokes equations, 
the entropy four-current is decomposed as, S'^ — su'* -I- 
(3q^, where q^ is the heat flux, s is the local rest frame 
entropy density and (3 = 1/T is the inverse temperature. 
These last two scalar quantities satisfy the fundamental 
relation of thermodynamics, s = (3{e+p), for matter with 
no conserved charge. Hence, in Ist-order theories the 
only way to satisfy the second law of thermodynamics, 
using a linear relationship between the thermodynamic 
force and flux, is to choose, 



T^NS 

Hjvs 



VO, 



n 



3 



(27) 
(28) 



where rj and ^ are positive coefhcients of shear and bulk 
viscosity respectively, while = df^u^^ is the expansion 
scalar. The Ist-order theories (contrary to 2nd-order the- 
ories) are known to have intrinsic problems attributed to 
the immediate appearance and disappearance of the ther- 
modynamic flux once the thermodynamic force is turned 
on or off. As shown by Hiscock and Lindblom [33| the 
linearized version of these equations propagate pertur- 
bations acausally, and even though initially these might 
be weak signals they may grow unbounded bringing the 
system out of stable equilibrium. 

To remedy some of these problems the 2nd-ordcr the- 
ory of Israel and Stewart was constructed Q , similarly to 
the non-relativistic theory by Miiller [34| . This was built 
around the assumption that the entropy four-current con- 
tains second order corrections in dissipation due to vis- 
cosity (here we disregard heat conductivity and cross 
couphngs) such that, 5'' = sm'' -|- Pq^" - {(3/2){(3o]l^ + 
3/327r^/2)u^, where (3q and (32 are thermodynamic coeffi- 
cients related to the relaxation times. Applying the law 
of positive entropy production and some algebra leads to 
the transport equations for the shear and bulk pressure. 



1 






(29) 
(30) 



where the relaxation time of bulk viscosity and shear are, 
Ttt = 277/^2 and Tn = C/^o- The above equations are re- 
ferred as the truncated Israel-Stewart equations, since 
terms involving the divergence of the flow field and ther- 
modynamic coefficients have been neglected compared to 
the equations by Israel and Stewart [8|. However, the cur- 
rent form of the transport equations already captures the 
essential features of relaxation phenomena which make 
the theory causal and stable. 

Another crucial difference between first and second or- 
der theories is in the mathematical structure of the equa- 
tions. In Ist-order theories the viscous corrections ap- 
pear linearly proportional to the divergence of the flow 



field, therefore they are more sensitive to fluctuations 
and the inhomogeneities in the flow field, see section 
IIVI In 2nd-order theory not only the coefficients of 
viscosity but also the thermodynamic coefficients need 
to be specified. For example the later parameters are 
know for a relativistic Boltzmann gas of massive parti- 
cles (Jj, /3o = 216(to/3)~'*/p, leading to a relaxation time 
of, Tn = 216(m/3)~^(C/p), which is divergent for a fiuid of 
massless particles, while the bulk viscosity coefficient is 
zero in that limit. The thermodynamic coefficient for the 
shear viscosity within the same substance is, (32 = 3/(4p), 
leading to the relaxation time of, r^ ~ (3/2) (77 /p). 

In passing we point out a few important facts regard- 
ing the bulk viscosity and its source. For a long time 
in the classical non-relativistic Navier-Stokes theory the 
purpose of bulk viscosity was controversial [35|. Even, 
Eckart in his pioneering work [3| was concerned with flu- 
ids without bulk viscosity. Israel was the first to show 
that the bulk viscosity of relativistic matter may not be 
unimportant [36|. This turned the attention mainly in 
cosmology, rendering bulk viscosity as the only possible 
form of dissipative phenomena, see for example 37, 3J[ 
and references therein for a thorough introduction. The 
bulk viscous effects are important in mixtures when the 
difference in property between the components becomes 
substantial. This might be due to the difference in cool- 
ing rates within the same type of substance or in a mix- 
ture between massive and effectively massless particles 
during a phase-transition [331 , (for other possible sources 
for bulk viscosity see [40|). Therefore in these situations 
bulk viscosity is used to describe a mixture effectively as 
a single fluid with a non-vanishing bulk viscosity coeffi- 
cient and relaxation time. 

It is also important to phenomenologically understand 
the relation between viscosity and relaxation time [3]. 
For example, bulk or volume viscosity appears when the 
system undergoes an isotropic expansion or contraction. 
If this happens at a relatively fast rate such that the sys- 
tem is unable to follow the change in volume and restore 
equilibrium in a short time, means that the relaxation 
time of the viscous pressure is long. On the opposite, 
if the system equilibrates almost immediately, than the 
corresponding relaxation time must be short. Hence it is 
also intuitive that large deviations from equilibrium can 
only be the consequence of large viscosity, while small 
departures from equilibrium result from small viscosity, 
assuming in both cases that the expansion rate is consid- 
erably small. It is also fundamental that the relaxation 
time must be shorter than the inverse of the expansion 
rate of the system, Ttt^h <C 1/9, otherwise the system 
will never be able to equilibrate and the fluid dynamical 
approach is unsuitable. 

In case of a relativistic Boltzmann gas, recalling the 
viscosity coefficients from eqs. (|27I28I) . we find that the 
relaxation times are 



T^NS 



^eq 



(31) 



Tn = 



in fTlNs\ 
\Pe, J ' 



(32) 



where d^ = 9/8 and we can only assume similarly to 
the relaxation time for the shear, that rfn is a dimen- 
sionless positive number on the order of unity. This also 
follows from the fact that both viscosities give birth to 
a local dissipative pressure, which for small dissipative 
corrections relax to the Navier-Stokes values. In case 
the dissipative pressure is comparable to the equilibrium 
pressure, the relaxation times become longer than the 
mean free time between collisions, thus the fluid dynam- 
ical approach may no longer be appropriate. 



III. THE NUMERICAL SCHEME 

Here we briefly review the basic principles of the un- 
derlying numerical scheme used in this work. The ex- 
plicit finite difference scheme called sharp and smooth 
transport algorithm (SHASTA) |41j] is a version of the 
flux corrected transport (FCT) algorithm. Detailed tests, 
simulations and comparisons to semi-analytical solutions 
have been performed with this algorithm in various situa- 
tions; in non-relativistic and relativistic perfect fluid dy- 
namics and magnetohydrodynamics; in the last decades 
p3 . [43 . m, |4^, \^ achieving confidence and wide usage. 

Before discussing the version of the algorithm in detail, 
let us rewrite the conservation equations (|16ll7ll8p and 
transport equations (|29I30|) in conservation form which 
makes it possible to treat all equations with the same 
numerical scheme. Due to similarity in form and effect in 
case of 1+1 dimensional expansion scenarios, we include 
only one type of viscosity and relaxation equation and 
refer to it as bulk viscosity in the following. Hence, 

dtR + d,iv,R) = 0, (33) 

dtE + d,{v.E) = -d,{v,P,), (34) 

dtM,+d,iv,M,) = -9,F,, (35) 

where R^ N^, E = r°o, M^ = r°^ and P^ = p + n or 
Pz = p + 11- Introducing a common notation, $ = 7$, 
for the auxiliary variables, tt = ^n and/or H = 7H, the 
relaxation equations (|29l30p can be rewritten^ in a form 
similar to the conservation equations. 






6*$. 



(36) 



where $, ^ns Sind t$ commonly denotes the shear pres- 
sure and its relaxation time and/or the bulk viscous pres- 
sure and its relaxation time. 

The above conservation and transport equations are of 
conservation type, written generally as. 



dtU + d,{vU)^S{t,z), 



(37) 



^ Another possibility would be to rewrite the relaxation equation 
as, dt^ + d4v,<S>) = -J- {<S>Ns -<!>) + {d,v,)'S>. 



where U = U{t, z) is one of the conserved quantities, 
V — v^is the velocity, and S{t, z) is the source term. The 
discretised conservative variable defined as an average in 
cell, j, at coordinate point, Zj, at discrete time level, t", 
is denoted by C/". Some of the source terms in our ex- 
amples contain differential operations, therefore they are 
represented as finite (second order) central differences, 
i.e., for spatial derivatives /S.Sj = (•S'JYi — S^_i)/2Az. 

In the SHASTA algorithm the velocity, the local rest 
frame variables and source terms, are computed and up- 
dated at half time steps, i.e., in At/2 time intervals. This 
requirement ensures second order accuracy in both space 
and time. In contrast, the conservative variables, U, used 
to advance the solution from time level n to n + I, are 
updated only once at the end of full time steps. In a 
given cell, J, this can be summarized formally as, 

f^«+l/2_^n(^n^^„^^n) ^ (38) 

[/"+! ^ U" fu", f"+l/2, S'"+l/2^ (39) 

In case of the relaxation equation, the source terms 
contain dynamical information on the divergence of the 
flow field in both space and time. Second order accu- 
racy in time can only be calculated in Af time inter- 
vals (if we use the time-split method), at time levels, 
n — 1/2, n,n+ 1/2, n + 1, . . ., where the time derivatives 
are, AS'" - (5" - S'^'^)/At, AS'"+i/2 - (S'"+i/2 _ 
S^~^^'^)/At, etc. This ensures better accuracy, (how- 
ever, the difference is rather small), than calculating the 
time derivatives as well as the source terms at full time 
steps only, i.e., only between time levels n and n + 1. 

The difference of primary variables in adjacent cells 
is denoted by Aj = C/JYi ~ ^j" °^ later also by A^- = 
Uj+i - ilj. The explicit SHASTA method [4l| at half- 
step as well as at full-step, first computes the so-called 
transported and diffused quantities. 



c/, 



(Q+ 



Q'-A,_i) 



(40) 



g_)c/;' + AiA5, 



where 



Q± = 



1/2 T e, 



1 ± (Cjil - €j) ' 
\ n+l/2 



(41) 
(42) 



and the Courant number is the ratio of time-step to cell- 
size, A = At/Az. A general requirement for any finite 
difference algorithm is to fulfill the so-called Courant- 
Friedrichs-Lewy (CLP) criterion, i.e., A < 1, related 
to the stability of hyperbolic equations, otherwise the 
numerical solution becomes unconditionally unstable. 
Physically this expresses that, matter must be causally 
propagated at most Az — At distance into vacuum. Por 
SHASTA, A < 1/2, while in this paper we use a smaller 
value, A = 0.4. Here we note that since the numerical al- 
gorithms average the transported quantities over a cell. 



part of the matter is acausally propagated over (1 — A)Az. 
This is a purely numerical artifact called prediffusion. 

The time-advanced quantities are calculated removing 
the numerical diffusion by subtracting the so-called an- 
tidiffusion fluxes, A, from the transported and diffused 
quantities such that, 



u: 



n+l 



U, - A, 



', -J -J . ^j-i • (43) 

Here we have defined the flux corrected antidiffusion flux 



A, 



cr," mm 



0, max 



\^2~^2 



+i,|A,|,a,A,_i ,(44) 



where the 'phoenical' antidiffusion flux^ is. 



A, = ^ 
^ 8 

CTj = sgn(Aj 



^j-^(A,+i-2A,+A,_i) 



, (45) 
(46) 



The so called mask, Aad-, is introduced to regulate the 
amount of anti-diffusion [47|- The algorithm tends to 
produce small wiggles, due to the fact that in the an- 
tidiffusion step one removes to much diffusion, therefore 
adjusting the mask one can suppress this artifact lead- 
ing to a more stable and smoother solutions. However 
the drawback is that, by reducing the antidiffusion we 
increase the numerical diffusion causing larger prediffu- 
sion and entropy production even in perfect fluids! This 
step is unavoidable in numerical algorithms where due to 
discretization the differential equations are truncated al- 
ready at leading order and without additional but purely 
numerical corrections lead to unstable solutions. Within 
the numerical framework this is called numerical dissipa- 
tion, or since it acts similarly to the physical viscosity it 
is also called artificial or numerical viscosity [481] . 

In case we want to model physical viscosity one has to 
keep in mind that the there is already a small numerical 
viscosity in the algorithm, which has to be estimated (for 
example by measuring the entropy production in case of 
a perfect fluid) and taken into account. This leads to 
a total effective viscosity which is larger than the one 
we explicitly include. Obviously, things are not as sim- 
ple, since the numerical viscosity and numerical diffusion 
contain linear and non- linear parts [46], and its effect 
strongly depends on the grid size, initial condition and 
flux limiters we use. Therefore it is a question of numer- 
ical analysis and extensive testing to reveal the effect of 
numerical viscosity. Some can be found in the original or 
related publications of the numerical schemes. 

It is important to remember that SHASTA is a low im- 
plicit viscosity algorithm, and conserves energy (and mo- 
mentum) up to 5-digits, but produces entropy in the case 
of a perfect fluid roughly between 0.5% — 5%, depending 
on the initial setup, antidiffuison flux, mask coefficient 



and physical situation. The lower value was found in the 
studies we are going to show in the next section using 
a mask of Aad = 0.8, while the error is less than 0.2% 
using the standard value, Aad = 1, after 200 time-steps. 
The large entropy production was found in the case of a 
3D grid with the same proportions, cell size, number of 
time-steps and reduced mask coefficient in case of 1-f 3 di- 
mensional expansion into vacuum of an initially constant 
energy sphere. 

To determine the bulk viscous pressure one also has 
to calculate the expansion scalar, 6{t,z). One possibil- 
ity is to take the standard form, used in this work, us- 
ing a second order accurate central difference formula, 
= dt-f + dziv^j) = j^ivzdtVz + d^v^). The other 
form can be expressed from the conservation of energy, 
Uyd^T^" = 0, or in case we also have conserved charge, 
from the continuity equation, 9^iV^ — 0, leading to, 9 = 
-^{dte + v^d^e) /{e+p + 11 + 71) = -j {dtn + Vzdzn) /n. 

The numerical differentiation of the velocity field in- 
troduces obvious numerical problems which we need to 
address. In particular, finite differences of the velocity 
field in adjacent cells, usually fluctuate due to numerical 
noise. Since we solve a set of non-linear coupled par- 
tial differential equations, these may become uncontrol- 
lable. To make the numerical expansion rate smoother 
we found that an additional five-point stencil smoothing 
is necessary^. However, the maximum number of neigh- 
boring cells to include is restricted by the Courant num- 
ber, otherwise we acausally propagate information into 
the neighboring cells. In our example the maximal num- 
ber of these cells are two to the right and two to the left, 
hence the five-point stencil. 

In the 1st- and 2nd-order theories the dissipative pres- 
sure, |H| < Peg, must be smaller than the equilibrium 
pressure. If the correction to the equilibrium pressure is 
small, the system will continue to expand with a lower 
effective pressure but the overall behavior should not 
change considerably from the perfect fluid limit. How- 
ever, at different parts of the system the local expansion 
rate may become very large (for example, in the tran- 
sition region to vacuum) and generate large dissipative 
corrections. This threshold is given by the equilibrium 
pressure locally, at least in the Ist-order theory. Hence, 
even though the physical situation may encounter larger 
values of the bulk pressure, we choose to keep this max- 
imum. The upper bound imposed on the bulk pressure 
leads to an upper bound for the local expansion rate, 
9max — Peq/C- I^. othcr words, the Navier-Stokes bulk 
pressure is defined to be U^s = —C,9 iov 9 < 9 max and 
j-jmaa; _ —(,9 max for 9 > 9 max- In the latter region the to- 
tal pressure vanishes and the acceleration therefore stalls. 

We keep the above convention also for 2nd-order the- 



The explicit antidiffusion flux |4ll |. Aj = AadAj/S, leads to 
somewhat smoother results. 



^ In a loose sense this coarse-graining of the expansion rate may be 
viewed as providing a "mass" to fluctuations with wave-length 
on the order of the grid spacing. 



ory so that for very short relaxation times we exactly 
approach the Navier-Stokes limit. It is also important to 
mention that since we solve the relaxation equation in 
conservation form, the expansion rate appears explicitly 
in the truncated equations as well, similarly as in the full 
Israel- Stewart equations. Hence it is obviously necessary 
to find an upper bound for terms containing the expan- 
sion rate otherwise those terms may grow unbounded and 
destabilize the solutions. 



IV. RESULTS AND DISCUSSION 

In all numerical calculations, we have fixed the follow- 
ing parameters. The Courant number, A = 0.4, and the 
cell size is dz — 0.2 fm, therefore dt — 0.08 fm/c. The 
grid contains 240 cells, while the total number of time 
steps is, n — 200, which corresponds to At — 16 fm/c 
expansion time. The amount of anti-diffusion is reduced 
by 20%, i.e., Aad — 0.8, which leads to some prediffu- 
sion but returns smoother profiles. The thermodynamic 
quantities are given by the Stefan-Boltzmann relations 
where the degeneracy of massless particles is 5 = 16. In 
case of dissipative fluids, the bulk viscosity coefficient to 
entropy density is constant with values corresponding to 
small, (/s — 0.2, and to large, C/s = 1, ratios. 

In most situations the initial expansion rate is un- 
known, therefore the dissipative corrections are neglected 
at start. This may only last for one time-step, since after 
that the time-derivatives can already be calculated and 
dissipative corrections added. In particular cases such as 
the Bjorken scaling solution, the expansion rate can be 
inferred from the geometry, therefore it does not pose a 
problem. We will return to this issue in Sect. IIVEI 

The relaxation time is given similarly to eq. (|32p . 

therefore rn = — i ( tt^ 1 • We have checked the asymp- 

t' y Peg J 

totic limits of the relaxation equations. In case the relax- 
ation time is small, rn ~ dt, the effect of bulk viscosity 
is immediately felt by the system, therefore the system 
behaves as in case of Ist-order theories. For very large 
relaxation times, i.e., larger than the lifetime of the sys- 
tem, and small initial dissipation, the effect of viscosity 
is exponentially suppressed and the solution approaches 
the perfect fluid limit. 



A. Expansion into vacuum of perfect fluid 

In case of perfect fluids one of the analytical solu- 
tions which relates the thermodynamic properties of mat- 
ter and the type of fluid dynamical solution is the H-1 
dimensional expansion into vacuum. This is a special 
case of the relativistic Riemann problem describing one- 
dimensional time-dependent flow. The initial conditions 
are such that initially at i = half of the space, z < 0, 
is filled uniformly with fiuid at rest, ^(z, 0) = 0, with 
energy density, e(z, 0) = eg, while the positive half z > 



is (empty) filled with vacuum. 

One can show that for thermodynamically normal mat- 
ter, e.g., a massless ideal gas with the EOS, p(e) = c^e, 
where c^ = 1/3 is the speed of sound squared, the stable 
solution to the fluid dynamical equations is a simple rar- 
efaction wave. The rarefaction wave is a wave where the 
energy density decreases in the direction of propagation, 
but the profile of the fiow does not change with time as 
a function of the similarity variable. 



u(e) 



1 - v(e)Cs 



(47) 



Here we recall the analytic results for a perfect fluid in the 
forward light-cone [43| . The energy density as a function 
of the similarity variable for, —1 < ^ < — c^, is constant. 



e(0 = c/po, 



(48) 



while in the region — c^ < ^ < 1 the matter starts to 
rarefy and the energy density decreases as 



e(0 = c/po 



1 + cs 1 + C 



(l+c^)/(2c,) 



(49) 



The temperature can be inferred from standard thermo- 
dynamical relations and the EOS, leading to 



T(£) = To ('^ 
eo 



c^/(i+=^) 



(50) 



We also compare how well the fluid flow is reproduced by 
the numerical calculation, where the analytical solution 
as a function of energy density is. 



w(e) = 



1 - (e/eo)^^°/(i+^') 
1 -K (e/eo)2==/(i+^') 



(51) 



The velocity can be given as a function of the similarity 
variable as well from eq. (|T7)) . this is plotted in Fig. [I}l. 
The results for the expansion scalar calculated numeri- 
cally are shown in Fig. [TId. 

Fig. [2^ shows the temperature normalized by the ini- 
tial temperature, from eq. ((50| . while Fig. [2)d, shows the 
lab frame energy density normalized by the initial pres- 
sure as a function of the similarity variable. The thick 
line shows the exact solution for a perfect fluid (ES-PF), 
while the numerical solutions also for a perfect fluid are 
at t = 4 fm/c with thin dotted line, at t = 8 fm/c with 
thin dashed line, and at i = 16 fm/c with fifll thin line. 

Both figures compare the analytical solutions to the 
numerical solutions, to pinpoint how well the underlying 
numerical method reproduces the ES-PF. We see that 
the numerical solutions asymptotically approach the ex- 
act result, while also reduce the numerical prediffusion 
into vacuum. This is due to the fact that the initially 
sharp discontinuity smears out as the rarefaction wave 
covers an increasing number of cells. Due to prediffusion 
and coarse-graining, the expansion rate it is not zero for 
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FIG. 1: The exact solution in case of a perfect fluid (thick line) 
and numerical solutions (thin dotted line at f = 4 fm/c, thin 
dashed line at, f = 8 fm/c, and thin line at t = 16 fm/c) as 
function of the similarity variable ^(z, t) . (a) the flow velocity 
v(^); (b) the expansion scalar Q(X)- 
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FIG. 2: (a) the temperature normalized to the initial tem- 
perature, r(f)/ro; (b) the laboratory frame energy density 
normalized to the initial pressure r'"'(^)/po; both as function 
of the similarity variable ^(z, i). The paramteres are the same 
as in Fig. [l] 



^ > 1 as can be seen in both figures; this problem mostly 
affects the acausally propagated low density matter. For 



a perfect fluid the numerical results are smooth and very 
nicely reproduce (especially at later times) the exact re- 
sults. The larger deviations from the ES-PF around the 
boundary to vacuum is due to a somewhat large prediffu- 
sion caused by the reduced antidiffusion, while the devia- 
tions around — c^ are due to a larger numerical diffusion. 
A more thorough study of the expansion into vacuum of 
a perfect fluid can be found in Ref. [43| . 



B. Expansion into vacuum with small and large 
dissipation 

Here we analyze and study the behavior of dissipative 
fluids corresponding to the Ist-order and 2nd-order the- 
ories, and plot the numerical results next to the exact 
solution in case of a perfect fluid. In all figures, unless 
stated otherwise, all quantities are plotted as a function 
of the similarity variable; the ES-PF is plotted with dot- 
ted line, while the numerical solutions at t == 4 fm/c with 
thin dashed line and at t = 16 fm/c with continuous line. 
The upper bound for the bulk pressure is, n™°^ = peg, 

and the relaxation time is, rn = — ^ ( ^^^^^ J fm/c. The ra- 
tio of viscosity over entropy density corresponds to small 
C/s = 0.2 or large C/s = 1.0 values. 
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FIG. 3: The velocity profile for a perfect fluid ES-PF com- 
pared to dissipative fluids with C,l s = 0.2 (thin) and C,l s — 1.0 
(thick), as numerically calculated and plotted at i = 4 fm/c 
(dashed) and t = 16 fm/c (full); (a) Ist-order theory; (b) 
2nd-order theory. 

Fig. [3^ shows the velocity profile calculated from Ist- 
order theory, while Fig. [3}d from 2nd-order theory. We 
see that for Ist-order theory the velocity is fluctuating 
with an increasing frequency in time. Initially the am- 



plitude and the wavelength of the fluctuations is large, 
however as the expansion proceeds these large fluctu- 
ations are damped and become smaller amplitude and 
smaller wavelength oscillations. This is partially due to 
that, that the initial discontinuity smears out in time 
and the problem is resolved on a larger grid. However, 
the damping is also due to the non-linear antidiffusion 
term (I44p in the numerical algorithm. In other words this 
is numerical viscosity, which always acts to smooth the 
fluxes, contrary to the numerical dispersion which acts 
exactly the opposite and produces ripples in the results. 
It should also be obvious that working on a larger grid 
with the same cell size would not improve the results! 

For bigger viscosity the numerical fluctuations become 
larger which is clearly a sign that the method fails and 
the numerical errors become uncontrollable. The results 
for 2nd-order theory are much smoother and show that 
there is a relevant change in the velocity profile, which is 
an outcome of the large dissipative pressure. This inter- 
esting phenomena appears near the edge of the matter 
where due to a large bulk pressure contribution the ef- 
fective pressure essentially vanishes, therefore that part 
of the system stops to accelerate but due to inertia it 
keeps moving forward with constant speed, hence form- 
ing a constant velocity plateau. Note that, the small 
wiggle in the velocity profile for large dissipation, visible 
for example on Fig. [Sb, is an artifact of the phoenical 
antidiffusion fiux. 
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FIG. 4: The calculation of the expansion scalar at different 
time-steps in case of a dissipative fluid with (/s = 0.2. (a) 
Ist-order theory; (b) 2nd-order theory. 

Figs. S^ andUjD shows the expansion scalar calculated 
numerically at different time-steps for both theories. We 
can see that in case of Ist-order theory, the expansion 




0.5 



FIG. 5: The effective pressure normalized by the initial pres- 
sure with dissipation proportional to (/s — 1. (a) Ist-order 
theory; (b) 2nd-order theory. 



scalar is highly fluctuating, while in 2nd-order theory it is 
much smoother. Both calculations reflect the space-time 
inhomogeneity of the flow field amplified by the velocity. 
For a better understanding we have plotted the effec- 
tive pressure as a function of the similarity variable in 
Fig. [5l for both theories. One can see that for early 
times the expansion rate and therefore the bulk viscous 
pressure is largest. When the effective pressure drops (to 
zero), the acceleration of matter is reduced (the matter 
flows with constant velocity). However, the expansion 
rate will also decrease later, reducing the viscous pres- 
sure and therefore the velocity will continue to increase. 
The important thing to remember is that the speed of 
sound decreases due to dissipation and the effective rar- 



efaction speed [49| can be given as, c; 



eff 



(p + n)/e. 



It is also interesting to remark that even though the 
bulk viscosity is large at places, the effective pressure 
(and therefore the equilibrium pressure) may be larger 
than in the case of perfect fluid, due to the fact the ther- 
modynamic state of the system is influenced by entropy 
production and slower cooling. This is why the lab frame 
energy density decreases slower, see Fig. [71 

Fig. [S] shows the temperature profile plotted at differ- 
ent time-steps for both theories. As in the case of velocity 
the presence of viscosity is observable in the overall re- 
duced cooling of matter. The increase in the expansion 
rate increases the dissipation which in turn slows the ex- 
pansion, thus the matter cools at a slower rate. On the 
other hand in Ist-order theories, the visibility of this ef- 
fect is much more reduced due to numerical problems. 

The laboratory frame energy density normalized to the 
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FIG. 6: The temperature normalized by the initial tempera- 
ture, (a) Ist-order theory; (b) 2nd-order theory. The param- 
eters are the same as in Fig. O 
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FIG. 7: The laboratory frame energy density normalized to 
the initial pressure, (a) Ist-order theory; (b) 2nd-order the- 
ory. The parameters are the same as in Fig. |3l 



initial pressure is presented in Fig. [T] Since this quantity 
is proportional to the fourth root of the temperature di- 
vided by the initial temperature, it is much less affected 
by the fluctuations and prediffusion in both cases. Based 
on the previous arguments the deviation from the ES-PF 



are noticeable for larger values of the similarity variable, 
where the effect of dissipation is most pronounced. This 
plot also confirms that the effective pressure drops to 
zero, however since the matter has a finite energy density, 
temperature, and velocity, the laboratory energy density 
is not zero at those places. As soon as the expansion 
rate decreases, the finite albeit small effective pressure 
will continue to expand the matter into vacuum. Fur- 
ther comparisons between the Ist-order and 2nd-order 
theories can be found in the Appendix. 



C. Expansion into vacuum with a soft EOS 

To further investigate the behavior of matter with 
large viscosity we have used a relatively soft EOS, where 
the speed of sound squared is c^ = I/I5, while keep- 
ing the other parameters intact. Using a soft equation 
of state reduces the pressure and the pressure gradi- 
ents in the system, while the relaxation time increases, 
Tn ~ (1 -I- c^)/(Cs T), inversely with local speed of sound 
and temperature. 




FIG. 8: The velocity profile for a perfect fluid, (ES-PF), com- 
pared to a dissipative fluid with C^/s = I, using a soft EOS. 
(a) fst-order theory; (b) 2nd-order theory. 

Fig. [5^ shows the velocity profile for Ist-order theory 
while Fig. [8)3 for 2nd-order theory, both in case of a soft 
EOS. In comparison to Fig. [31 the velocity profiles are 
much improved and the constant velocity part of velocity 
profile is clearly visible even from the Ist-order theory. 
It is clear that the algorithm works much better, pro- 
ducing overall smooth results. We have also checked our 
algorithm with a hard EOS, e = p, which proved that 
overshoots and oscillations become enhanced compared 
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FIG. 9: The effective pressure normaiized by ttie initiai pres- 
sure for a soft EOS. (a) Ist-order tfieory; (b) 2nd-order theory. 



to the standard case (c^ = 1/3) and the results became 
less reliable. 

Fig. [9l shows the effective pressure normalized by the 
initial pressure, similar to Fig. [5l but with a soft EOS. 
Once again the results are good, especially in the case of 
2nd-order theory. This is due to fact that a softer EOS 
not only reduces the thermodynamic pressure compared 
to the standard EOS but also decreases the dissipative 
pressure. It is also important to understand that in this 
case the dissipative effects act on an overall wider scale 
as function of ^, but even though the effect of dissipation 
is immediately added, the results are much smoother be- 
cause the dissipative pressure is also reduced. 



D. Expansion into vacuum of a matter with 
temperature dependent bulk viscosity 

A interesting and relevant question to study, is the ex- 
pansion of dissipative matter with a temperature depen- 
dent bulk viscosity. To model this property, we assumed 
that bulk viscosity acts in the close vicinity of a specific 
or critical temperature, T = Tc(l ± 0.02), otherwise it is 
zero. Thus the bulk viscosity coefficient is. 



C = Coe(r - i.02r,) (i - e(r - o.gsr,)) 



(52) 



Here the choice of critical temperature is, Tc — 2To/3, 
where Tq is the initial temperature, and ^o = s is the bulk 
viscosity coefficient. The calculations are done for 2nd 
order theory including the above temperature dependent 
bulk viscosity. 



FIG. 10: The velocity profile (a), and the numerical calcula- 
tion of the expansion scalar (b) , using a temperature depen- 
dent bulk viscosity for 2nd-order theory. 
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FIG. 11: The parameters are the same as in Fig. 1101 (a) the 
temperature normalized to the initial temperature, T{(,)/To; 
(b) the laboratory frame energy density normalized to the 
initial pressure r'"'(5)/po. 



Figs. [TOk and [TOb shows the velocity of the matter 
and the expansion rate as a function of the similarity 
variable. For the same setup, in Figs. [TTb andfTTb the 
temperature normalized by the initial temperature and 
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the laboratory frame energy density normalized by the 
initial pressure is shown. 

When the temperature falls in the respective regime, 
the viscosity over entropy ratio rises suddenly, to (q/s — 
1, while otherwise the dissipation is switched off. In our 
case this manifests itself as almost constant velocity tem- 
perature, pressure, and energy density plateau, located 
roughly between < ^ < 0.2. Because the effective 
pressure decreased suddenly, the matter is slowed down 
considerably, until the matter cools below the predefined 
temperature (although much slower), then the system 
will suddenly find local thermal equilibrium and continue 
to accelerate. This is apparent in all plots. This type of 
studies may be relevant in case of a phase-transitions 
where the temperature dependent viscosity modeling be- 
comes necessary |5G |. 



E. The Bjorken solutions for perfect and 
dissipative fluids 

In this section we test how well the numerical calcula- 
tions reproduce one-dimensional dissipative scaling flow. 
The relaxation time was kept constant, rn = 1 fm/c, 
which does not effect the final outcome qualitatively. 

We first recall the 1-1-1 dimensional Bjorken scaling 
solution for perfect [5l| and dissipative fluids 5]. The 
equations follow from the conservation law d^/T^^'^ — 0, 
and the second law of thermodynamics, under the as- 
sumption that the matter expands longitudinally with a 
flow velocity v = z/t in a boost invariant manner. To 
transform the partial differential equations into simple 
differential equations (using the assumption of boost in- 
variance), one carries out a coordinate transformation 
from (i, z) to (t, 77), where r = Vi^ — z^ is the proper 
time and il = \ log [{t + z)/{t — z)] is the space-time ra- 
pidity. Therefore, the truncated Israel-Stewart equations 
for the energy and bulk pressure are, 



de 
d7 



1 



{e+p + n) 



(n 



NS 



m 



(53) 
(54) 



where Uns — ^Cl'''-, a-nd the effective pressure satis- 
fies, d{p + n)/d77 = 0. The equations of perfect fluid 
dynamics are obtained when the dissipative bulk pres- 
sure is zero, n(T) — 0. The relativistic Navier-Stokes 
equation is given by eq. (|53p alone, since bulk pressure 
n(T) = nAr5(r) is given algebraically. 

In the Bjorken picture the system is infinitely elon- 
gated in rapidity. Since our SHASTA code is written in 
standard space-time coordinates (i, z), we have to deter- 
mine initial values on a i = io surface and the fluid must 
have a finite length (due to the finite grid), —Z[)<z<zq. 
This is done as follows. 

We first solve eqs. (I53l54p using a 4th-order Runge- 
Kutta solver for all times t > tq — 1 fm/c. Together 
with the assumption of boost invariance, this determines 



the hydrodynamic fields in the entire forward light-cone. 
We then repeat the solution with the SHASTA partial 
differential equation solver with initial conditions at io = 
zo-l- Az = 6 f m set by the Runge-Kutta solution. We note 
that in the SHASTA solution the system has a boundary 
since the velocity at zq is close to but slightly less than 
the velocity of light (zq is smaller than t^ by one grid 
spacing) . 




FIG. 12: The exact velocity vsj (full line) and expansion 
rate Qbj (dash dotted line), compared to the numerically 
calculated velocity v (dotted line) and expansion rate O (full 
line), for a dissipative fluid with C,/s = 0.2 after At = 8 
fm/c evolution, (a) Ist-order theory; (b) 2nd-order theory 
with rn = 1 fm/c constant relaxation time. In both cases the 
initial value for the bulk pressure is given by the Navier-Stokes 
value. 

In the Ist-order theory the initial value for the bulk 
pressure is Uns = — Co/''o, which can be limited by the 
equilibrium pressure, i.e., the dissipative pressure should 
not be larger than the initial equilibrium pressure oth- 
erwise the system becomes unstable. In the 2nd-order 
theory, we take the initial value of the bulk pressure to 
be the same as in the Ist-order theory. Using this ini- 
tial condition allows for a direct comparison since both 
theories start from the same initial values. Moreover, 
on physical grounds, if tq is to be interpreted as the on- 
set of hydrodynamic behavior (thermalization time) than 
the Reynolds number at tq should be stationary, i.e., it 
should neither grow nor decrease; this implies that the 
initial value for the bulk pressure should be close to that 
given by the Ist-order approach [53 |. 

We can see in Fig. [T2]that the flow velocity and expan- 
sion rate are fairly well reproduced for the 2nd-order the- 
ory while the results for the Ist-order theory show large 
numerical artifacts already at a few fm distance from the 
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FIG. 13: The effective pressure normalized by the initial pres- 
sure, (a) Ist-order theory; (b) 2nd-order theory. The param- 
eters are the same as in Fig. [T^ 



center. This is correlated with the coarse graining since 
the results improve on finer grids^. We can also observe 
the large predifFusion into vacuum due to reduced antid- 
iffusion similarly to the Riemann wave discussed above. 
The fluctuations in the velocity are also visible in the 
expansion rate and in the pressure shown in Fig. 1131 

We have also tested how well the algorithm solves the 
relaxation equation for the bulk pressure. This is impor- 
tant to know, since the SHASTA was specially designed 
to solve conservation equations. We have plotted the evo- 
lution of the bulk pressure in the central cell, which has 
a velocity w « 0, next to the results of a standard fourth- 
order Runge-Kutta solver. The comparison in case of 
Ist-order theory is given in Fig. [Hk . Here ^xi^rk) with 
full line is the Runge-Kutta solver, while the result of 
SHASTA is Hi, with dashed line. The curves match rea- 
sonably well, deviations are due to the overestimate of 
the local expansion rate. 

For 2nd-order theories our standard choice for the ini- 
tial value of the bulk pressure is given by the Navier- 
Stokes value; this corresponds to ^2{rk) ^-nd ^2 in 
Fig. [Tib . We have also compared to the case when the 
system starts from equilibrium^ (n(io) = 0), see the 



* In view of forthcoming applications to three-dimensional prob- 
lems we only consider grids with at most a few hundred grid 
cells. 

^ Note that, as mentioned above, for such type of initial conditions 
the dissipative corrections initially grow very rapidly, which is 
not a physically plausible scenario 1531 . 



FIG. 14: The evolution of the bulk pressure in a central cell. 
(a) Ist-order theory; (b) 2nd-order theory. The initial condi- 
tion for the bulk pressure is given by the Navier-Stokes value 
for. III, 112, ^\(RK) and ^2{rk)- The evolution of the initially 
equilibrated system are for, lis and ^^(rk)- 



curves $ 



Z{RK) 



and $•:(. In both cases the SHASTA result 



is very good. There are some small deviations, however, 
this is due to the difference between the calculated and 
analytical expansion rates. We have checked our calcu- 
lations using the exact expansion rate Q = 1/t and in 
this case not only the 2nd-order calculations but also the 
Ist-order ones are smooth and very accurately match the 
Runge-Kutta results. 



V. SUMMARY AND OUTLOOK 

In this work we have focused on testing numerical so- 
lutions of first and second order theories of relativistic 
hydrodynamics with bulk viscosity using the SHASTA 
flux corrected transport algorithm. This is a rather effi- 
cient and fast algorithm for solving causal fluid transport 
on a fixed grid; it provides accurate solutions of ideal hy- 
drodynamics with minimal numerical viscosity and pred- 
iffusion and can be easily adapted to multi-dimensional 
problems. In fact, the algorithm can also be used to 
solve the relaxation equations of the 2nd-order approach 
simultaneously with the conservation equations without 
resorting to other numerical schemes, which may reduce 
the computational time and complicate the problem and 
its implementation further. 

The Ist-order theory of viscous hydrodynamics pro- 
vides the proper description of long- wavelength, low- 
frequency density waves in a fluid. The 2nd-order theory 
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introduces relaxation equations for the dissipative fluxes 
thereby maintaining causality. Its solutions converge to 
those of the Ist-order theory over time scales larger than 
the relaxation times. It has been argued that these re- 
laxation times might be on the order of the microscopic 
time scales in the problem and that the 2nd-order the- 
ory is therefore no better approximation to the dynamics 
than the 1-st order approach. 

Our numerical solutions with the SHASTA algorithm, 
however, indicate that the accuracy and stability of the 
solutions of the 2nd-order theory is significantly better 
than in the Ist-order theory, even if the calculated local 
expansion rate is smoothed over a few fluid cells: the so- 
lution of the Riemann wave with viscosity in the Ist-order 
approach produces oscillations which are absent from the 
2nd-order theory. This observation holds for virtually 
any amount of dissipation. Also, the numerical problems 
encountered in the Ist-order approach get milder if the 
speed of sound is smaller (which reduces the acceleration 
of the flow) but worse if the equation of state is stiff. 

These observations are valid in 2 or 3 dimensional cases 
[53|, therefore in conclusion, we believe that for general 
purpose codes the 2nd order theory is not only more gen- 
eral, but also more stable and reliable even numerically. 
Although using the 2nd-order theory it is computation- 
ally more intensive since the dissipative quantities have 
to be propagated in time, its implementation into some 
existing numerical codes which solve hyperbolic partial 
differential equations in conservation form, does not re- 
quire more effort than adding the Ist-order corrections. 

Regardless on the difference between numerical meth- 
ods, using the same initial conditions and correspond- 
ing physical quantities, all numerical results should be 
very closely the same. Unfortunately, in case of dissipa- 
tive fluid dynamics there are only a few simple solutions 
where the numerical accuracy can be tested in great de- 
tail. However, taking into account that actual applica- 
tions of relativistic fluid dynamics in modeling relativistic 
heavy-ion collisions need several other crucial approxi- 
mations, (introducing additional uncertainties and pa- 
rameters linked to the fluid dynamical calculation, such 
as initial conditions and freeze-out), it is of great impor- 
tance that the numerical fluid dynamical methods should 
be very carefully investigated, tested and documented in 
various situations. To our knowledge, in most publica- 
tions this topic is rather forgotten and/or undisclosed. 
The other reliable possibility and recommendation would 
be to check the fluid dynamical codes against kinetic the- 
ory [5J,[53, which on the other hand would also 'validate' 
transport codes in the fluid dynamical regime. 

Note added: During the preparation of this manuscript 
we became aware of the very recent work by the Brazilian 
group 56, 57], on the shock propagation and stability in 
causal 1+1 dimensional dissipative hydrodynamics, us- 
ing the smoothed particle hydrodynamics (SPH). This 
important problem was also investigated by us with the 
relaxation method presented in this paper, and lead to 
very promising agreement with kinetic theory [58| . 
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APPENDIX A 

Here we approximately extract and quantify the nu- 
merical errors and uncertainties as a function of time and 
viscosity. Since the exact solution for the Riemann prob- 
lem with dissipation is unknown, the deviation from an- 
alytic solutions cannot be measured. However, since the 
2nd order solutions converge to the Ist-order solutions 
at late times, we may quantify the error and deviations 
between the outcomes. 
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FIG. 15: The relative difference between the Ist-order and 
2nd-order theory as a function of time, measured by the above 
integral, for values of viscosity given in the figure. 

Therefore, we introduce the following relative devia- 
tion (point by point) between the flow velocities of matter 
and measure it by the following integral. 



dit) 
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(55) 



rl.5 



where the normalization factor, N — J_'^ d^\v2nd{0\- 
Note that, the Ist-order and 2nd-order theories start 
from different initial values, hence the initial deviation. 
We also observe that increasing the viscosity the results 
start to diverge at late times, signaling the instability and 
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large errors in the numerical solution of the Ist-order the- 
ory. This was apparent in all figures in Sect. IIVBI 
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FIG. 16: The velocity profiles as calculated numerically at 
t — 16 fm/c for a dissipative fluid with C^/s = 0.1 in 1-st order 
theory using different antidiffusion mask coefficients (ad) as 
shown on the figure. The initial condition in subplot (a) is 
V = 1 for 2 > 0, guarantees a continuous solution at the 
boundary to the vacuum. In (b) we have v — for z > 0. 

Next we discuss and show further examples and test 
results using different initial conditions and antidiffusion 
mask coefficients. Using the purely conventional initial 
value for the velocity of vacuum, v = 1 for z > 0, we 
make sure that the fluid dynamical solution is continuous 
at the boundary to vacuum. 

This specific test ensures that the oscillations in the 
1-st order theory are not caused by the sharp boundary 
and that the oscillations propagate outwards and not in- 
wards from the vacuum, see Fig. 116b . Therefore, this 
effect is mainly due to two things. First the higher order 
derivatives of the flow velocity appear in the equations 
and therefore even very small fluctuations in the flow field 
are enhanced and couple back into the solution. The sec- 
ond is a purely numerical problem which unfortunately 
affects of the numerical scheme and its accuracy since the 
SHASTA algorithm was not explicitly designed to solve 
the relativistic Navier-Stokes equations. 

We also see the effect of a further numerical artifact 
namely reducing the antidiffusion coefficient by 10%, 20% 
and 30% not only increases the entropy in the system, but 
also gives result to non-linear changes and differences in 
the solutions. The standard version of the algorithm uses 
a mask coefficient equal to 1. We can see in Figs. [TBI and 
[T7b that using a 20% smaller mask is reasonably close 
to the solution with the standard value of the mask, but 
more importantly the results become much smoother. 
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FIG. 17: The velocity profile for a Bjorken expansion in Ist- 
order theory for ("/s = 0.2 after Ai = 8 fm/c evolution. On (a) 
the un-smoothed vs. the smoothed vsm velocity profile with 
ad=0.8 mask coefficient. The continuous boundary condition 
for the velocity with various mask coefficients on (b). 



In Fig. [T7b we have plotted the velocity profile calcu- 
lated with a 20% reduced antidiffusion for a smoothed 
and un-smoothed expansion rate. We can see that the 
smoothing affects the solution positively leading to even 
less prediffusion into vacuum in this particular case. 



APPENDIX B 

Here we recall the solutions for the expansion into vac- 
uum in case of a perfect fluid following [591] . Introducing 
the similarity variable, ^ = z/t, the partial derivatives 
transform as, dt = — (C/0 {d/d^) and dz = (1/0 {d/d^). 
Therefore the equation for the energy and momentum in 
terms of the rest frame quantities becomes. 



dv 



(56) 



(57) 



+ 7^(6 + P) [{v - 0(1 + 2v^7^) + v] — ^0. 

Using the standard EOS, P = c^e, the vanishing de- 
terminant of the above system of equations leads to the 
expression for the characteristic variable, ^ — j^^- The 
correct solutions implies eq. (j47p . hence we are lead to 
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the following trivial equation, 

(1 + cl) dv 



de 
e 



(1-^2) ' 



(58) 



which with the corresponding initial conditions given in 
sect. IIV Al leads to the results presented before. 

Viscosity is introduced as, P — c^e + H, where the 
expansion rate is. 



5^u^ 






(59) 



Using the expression for the dissipative pressure in 1st- 



order theory, dP/d^, leads to terms containing, d^v/d^"^ ^ 
(dw/(i^)2, d(/d£, and dt'^/d^ = —t/^, therefore even in 
this simple case the exact solution is unknown. In 2nd- 
order theory the relaxation equation is. 



(v-Ol 



de, m d£, 



n, 



Tn 



(60) 



while the derivative of the pressure reduces to, dP/dS, = 
c^de/d^+dH/de,. In conclusion we see that for dissipative 
fluids the equations depend explicitly on the time in the 
local rest frame and the similarity of the flow is broken. 
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